import os
import sys
from cmath import cos
from lib2to3.pgen2.pgen import DFAState

import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
from scipy import interpolate


def GUM_RayD(
    data_L2,
    data_L0,
    param="D",
    z0=30,
    dens0=0.01841,
    u_dens0=0,
    tau=1e-9,
    u_tau=0,
    alt_angle=0,
    location_altitude=0,
    height_resolution=1,
    time_resolution=3600,
    pulse_repeat=30,
    k=1,
):
    # data_L0：L0级数据，第1列为高度，第2列为光子数
    # data_L2：L2级数据，第1列为高度，第2列为密度
    # tau：探测器渡越时间
    # u_tau：探测器渡越时间波动度
    # pulse_repeat：激光器重频
    # time_resolution：L2级积分时间
    # location_altitude：系统所在高度
    # z0：参考密度所在高度
    # dens0：参考密度
    # uncertain_densityRef：参考密度波动度
    # tau = 1e-9
    # u_tau = 1e-9 * 0.01
    # time_resolution = 1200
    # pulse_repeat = 30
    # z0 = 30  #参考密度所在高度
    # dens0 = 3.64686762151e+017  #参考密度
    # u_dens0 = 0  #参考密度波动度
    # location_altitude = 1.14
    z = data_L2[:, 0]
    dens = data_L2[:, 1]
    # 海拔高度转化和距离门内光子数累加
    data_L0_s = get_bins_data(
        data_L0, z, location_altitude, alt_angle, height_resolution
    )
    # 背景噪声获取
    B, B_mean_num = get_backgroud_noise(data_L0_s)
    B = B[0]
    # 探测范围截取
    Nr = get_target_data(data_L0_s, z)[:, 0]  # 目标高度光子数
    Nr0 = get_target_data(data_L0_s, z0)[0]  # 参考高度处光子数
    # z = z[:np.size(Nr)]
    # dens = dens[:np.size(Nr)]
    # 光子数校正过程中的灵敏系数
    Nc1 = Nr
    Nb1 = Nc1 - B
    allSize = time_resolution * pulse_repeat  # L2数据累加的脉冲数
    c_Nb1_tau1 = Nc1**2 / allSize
    c_Nb1_Nr = Nc1**2 / Nr**2
    c_Nb1_B = 1
    Nc2 = Nr0
    Nb2 = Nc2 - B
    c_Nb2_tau2 = Nc2**2 / allSize
    c_Nb2_Nr0 = Nc2**2 / Nr0**2
    c_Nb2_B = 1
    # 反演过程中灵敏系数
    c_dens_Nb1 = dens / Nb1
    c_dens_Nb2 = dens / Nb2
    c_dens_dens0 = dens / dens0
    # 光子噪声不确定度分量
    u_Nr = np.sqrt(Nr)
    u_Nr0 = np.sqrt(Nr0)
    c_Nr = c_dens_Nb1 * c_Nb1_Nr
    c_Nr0 = c_dens_Nb2 * c_Nb2_Nr0
    u1 = np.abs(u_Nr * c_Nr)
    u2 = np.abs(u_Nr0 * c_Nr0)
    # 背景噪声不确定度分量
    u_B1 = np.sqrt(B / B_mean_num)
    c_B1 = c_dens_Nb1 * c_Nb1_B
    c_B2 = c_dens_Nb2 * c_Nb2_B
    c_B = c_B1 + c_B2
    u3 = np.abs(u_B1 * c_B)
    # 渡越时间不确定度分量
    c_tau1 = c_dens_Nb1 * c_Nb1_tau1
    c_tau2 = c_dens_Nb2 * c_Nb2_tau2
    u4 = np.abs(u_tau * c_tau1)
    u5 = np.abs(u_tau * c_tau2)
    # 参考密度不确定度分量
    u6 = np.abs(dens0 * u_dens0 / 100 * c_dens_dens0)
    # 合成标准不确定度
    u = np.sqrt(u1**2 + u2**2 + u3**2 + u4**2 + u5**2 + u6**2)
    # 扩展不确定度
    U = k * u
    # 不确定度相对百分比
    U, u, u1, u2, u3, u4, u5, u6 = np.round(
        (U, u, u1, u2, u3, u4, u5, u6) / dens * 100, 3
    )
    result = np.vstack((z, dens, U, u, u1, u2, u3, u4, u5, u6)).T
    if param == "D":
        unit = "kg/m^3"
    else:
        unit = "1/cm^3"
    header = [
        "高度(km)",
        "密度({})".format(unit),
        "U(%)_GUM扩展不确定度",
        "u(%)_标准不确定度",
        "u1(%)_探测点光子噪声",
        "u2(%)_参考点光子噪声",
        "u3(%)_背景噪声",
        "u4(%)_探测点渡越时间",
        "u5(%)_参考点渡越时间",
        "u6(%)_参考密度",
    ]
    result_df = pd.DataFrame(result, columns=header)
    # result_df.drop('U(%)_GUM扩展不确定度', axis=1, inplace=True)
    return result_df


def GUM_RayT(
    data_L2,
    data_L0,
    z0=30,
    dens0=0.01841,
    u_dens0=0,
    z1=80,
    T1=209.6083,
    u_T1=0,
    tau=1e-9,
    u_tau=0,
    alt_angle=0,
    location_altitude=0,
    height_resolution=1,
    time_resolution=3600,
    pulse_repeat=30,
    k=1,
):
    # data_L0：L0级数据，第1列为高度，第2列为RayVHS
    # data_L2：L0级数据，第1列为高度，第2列为RayVHS
    # tau：探测器渡越时间
    # u_tau：探测器渡越时间波动度
    # pulse_repeat：激光器重频
    # time_resolution：L2级积分时间
    # location_altitude：系统所在高度
    # z0：参考密度所在高度
    # dens0：参考密度
    # uncertain_densityRef：参考密度波动度
    # tau = 1e-9
    # u_tau = 1e-9 * 0.01
    # time_resolution = 1200
    # pulse_repeat = 30
    # z0 = 30  #参考密度所在高度
    # dens0 = 3.64686762151e+017  #参考密度
    # u_dens0 = 0  #参考密度波动度
    # z1 = 80  #密度积分起始高度
    # T1 = 209.6083 #密度积分起始温度
    # u_T1 = 0  #参考温度波动度
    z = data_L2[:, 0]
    T = data_L2[:, 1]
    # todo
    # mask = z <= z1
    # z = z[mask]
    # T = T[mask]
    # 海拔高度转化和距离门内光子数累加
    data_L0_s = get_bins_data(
        data_L0, z, location_altitude, alt_angle, height_resolution
    )
    # 背景噪声获取
    B, B_mean_num = get_backgroud_noise(data_L0_s)
    B = B[0]
    # 探测范围截取
    Nr = get_target_data(data_L0_s, z)[:, 0]  # 目标高度光子数
    Nr0 = get_target_data(data_L0_s, z0)[0]  # 参考高度处光子数
    Nr1 = get_target_data(data_L0_s, z1)[0]  # 参考高度处光子数
    # np.savetxt('L0.txt', data_L0_s)
    # 光子数校正过程中的灵敏系数
    Nc = Nr
    Nb = Nc - B
    allSize = time_resolution * pulse_repeat  # L2数据累加的脉冲数
    c_Nb_tau1 = Nc**2 / allSize
    c_Nb_Nr = Nc**2 / Nr**2
    c_Nb_B = 1
    Nc0 = Nr0
    Nb0 = Nc0 - B
    c_Nb0_tau2 = Nc0**2 / allSize
    c_Nb0_Nr0 = Nc0**2 / Nr0**2
    c_Nb0_B = 1
    Nc1 = Nr1
    Nb1 = Nc1 - B
    # c_Nb1_tau2 = Nc1**2 / allSize
    # c_Nb1_Nr1 = Nc1**2 / Nr1**2
    # c_Nb1_B = 1
    # 反演过程中灵敏系数
    c_T_T1 = Nb1 / Nb * z1**2 / z**2
    c_T_Nb = -T / Nb
    c_T_Nb0 = T / Nb0
    # 光子噪声不确定度分量
    u_Nr = np.sqrt(Nr)
    u_Nr0 = np.sqrt(Nr0)
    c_Nr = c_T_Nb * c_Nb_Nr
    c_Nr0 = c_T_Nb0 * c_Nb0_Nr0
    u1 = np.abs(u_Nr * c_Nr)
    u2 = np.abs(u_Nr0 * c_Nr0)
    # 背景噪声不确定度分量
    u_B = np.sqrt(B / B_mean_num)
    c_B1 = c_T_Nb * c_Nb_B
    c_B2 = c_T_Nb0 * c_Nb0_B
    c_B = c_B1 + c_B2
    u3 = np.abs(u_B * c_B)
    # 渡越时间不确定度分量
    c_tau1 = c_T_Nb * c_Nb_tau1
    c_tau2 = c_T_Nb0 * c_Nb0_tau2
    u4 = np.abs(u_tau * c_tau1)
    u5 = np.abs(u_tau * c_tau2)
    # 参考温度不确定度分量
    u6 = np.abs(u_T1 * c_T_T1)
    # 合成标准不确定度
    u = np.sqrt(u1**2 + u2**2 + u3**2 + u4**2 + u5**2 + u6**2)
    # 扩展不确定度
    U = k * u
    # 不确定度
    U, u, u1, u2, u3, u4, u5, u6 = np.round((U, u, u1, u2, u3, u4, u5, u6), 4)
    result = np.vstack((z, T, U, u, u1, u2, u3, u4, u5, u6)).T
    header = [
        "高度(km)",
        "温度(K)",
        "U(K)_GUM扩展不确定度",
        "u(K)_标准不确定度",
        "u1(K)_探测点光子噪声",
        "u2(K)_参考点光子噪声",
        "u3(K)_背景噪声",
        "u4(K)_探测点渡越时间",
        "u5(K)_参考点渡越时间",
        "u6(K)_参考温度",
    ]
    result_df = pd.DataFrame(result, columns=header)
    # result_df.drop('U(K)_GUM扩展不确定度', axis=1, inplace=True)
    return result_df


def GUM_RayW(
    data_L2,
    data_L0,
    mode=0,
    c_V_RV=707,
    wavelen=532,
    u_wavelen=0.01,
    tau=1e-9,
    u_tau=0,
    alt_angle=0,
    location_altitude=0,
    height_resolution=1,
    time_resolution=3600,
    pulse_repeat=30,
    k=1,
):
    # c_V_RV：测风灵敏度
    # mode: 0单边缘，1双边缘

    # data_L0：L0级数据，第1列为高度，第2列为光子数
    # data_L2：L2级数据，第1列为高度，第2列为密度
    # tau：探测器渡越时间
    # u_tau：探测器渡越时间波动度
    # pulse_repeat：激光器重频
    # time_resolution：L2级积分时间
    # location_altitude：系统所在高度
    # z0：参考密度所在高度
    # dens0：参考密度
    # uncertain_densityRef：参考密度波动度
    # tau = 1e-9
    # u_tau = 1e-9 * 0.01
    # time_resolution = 1200
    # pulse_repeat = 30
    # z0 = 30  #参考密度所在高度
    # dens0 = 3.64686762151e+017  #参考密度
    # u_dens0 = 0  #参考密度波动度
    # location_altitude = 1.14
    z = data_L2[:, 0]
    V = data_L2[:, 1]
    # 海拔高度转化和距离门内光子数累加
    data_L0_s = get_bins_data(
        data_L0, z, location_altitude, alt_angle, height_resolution
    )
    # 背景噪声获取
    B, B_mean_num = get_backgroud_noise(data_L0_s)
    B = B[0]
    # 探测范围截取
    Nr = get_target_data(data_L0_s, z)  # 目标高度光子数
    Nr1 = Nr[:, 0]
    Nr2 = Nr[:, 1]
    # 光子数校正过程中的灵敏系数
    Nc1 = Nr1
    Nb1 = Nc1 - B
    allSize = time_resolution * pulse_repeat  # L2数据累加的脉冲数
    c_Nb1_tau1 = Nc1**2 / allSize
    c_Nb1_Nr1 = Nc1**2 / Nr1**2
    c_Nb1_B = 1
    Nc2 = Nr2
    Nb2 = Nc2 - B
    c_Nb2_tau2 = Nc2**2 / allSize
    c_Nb2_Nr2 = Nc2**2 / Nr2**2
    c_Nb2_B = 1
    # 反演过程中灵敏系数
    if mode == 0:
        # 针对RV = (N1-N2)/(N1+N2)的情況
        # if c_V_RV == None:
        #     RV = (Nb1 - Nb2) / (Nb1 + Nb2)
        #     c_V_RV = np.polyfit(RV, V, 1)[0]
        c_RV_Nb1 = 2 * Nb2 / (Nb1 + Nb2) ** 2
        c_RV_Nb2 = 2 * Nb1 / (Nb1 + Nb2) ** 2
    else:
        # 针对RV = N1/N2的情況
        # if c_V_RV == None:
        #     RV = Nb1 / Nb2
        #     c_V_RV = np.polyfit(RV, V, 1)[0]
        c_RV_Nb1 = 1 / Nb2
        c_RV_Nb2 = Nb1 / Nb2**2
    # plt.figure()
    # plt.scatter(RV, V)
    # plt.plot(RV, RV * c_V_RV + np.polyfit(RV, V, 1)[1])
    # plt.savefig('ResultPic/风速偏导拟合.png')
    # plt.close()
    c_V_Nb1 = c_V_RV * c_RV_Nb1
    c_V_Nb2 = c_V_RV * c_RV_Nb2
    # c_V_V0 = V / V0
    # 光子噪声不确定度分量
    u_Nr1 = np.sqrt(Nr1)
    u_Nr2 = np.sqrt(Nr2)
    c_Nr1 = c_V_Nb1 * c_Nb1_Nr1
    c_Nr2 = c_V_Nb2 * c_Nb2_Nr2
    u1 = np.abs(u_Nr1 * c_Nr1)
    u2 = np.abs(u_Nr2 * c_Nr2)
    # 背景噪声不确定度分量
    u_B1 = np.sqrt(B / B_mean_num)
    c_B1 = c_V_Nb1 * c_Nb1_B
    c_B2 = c_V_Nb2 * c_Nb2_B
    c_B = c_B1 + c_B2
    u3 = np.abs(u_B1 * c_B)
    u4 = u3
    # 渡越时间不确定度分量
    c_tau1 = c_V_Nb1 * c_Nb1_tau1
    c_tau2 = c_V_Nb2 * c_Nb2_tau2
    u5 = np.abs(u_tau * c_tau1)
    u6 = np.abs(u_tau * c_tau2)
    # 参考波长不确定度分量
    c_wavelen = 3e8 / (wavelen * 1e-9) / 2
    u7 = np.abs(u_wavelen * c_wavelen) * np.ones(np.shape(V))
    # 合成标准不确定度
    u = np.sqrt(u1**2 + u2**2 + u3**2 + u4**2 + u5**2 + u6**2)
    # 扩展不确定度
    U = k * u
    result = np.vstack((z, V, U, u, u1, u2, u3, u4, u5, u6, u7)).T
    header = [
        "高度(km)",
        "风速(m/s)",
        "U(m/s)_GUM扩展不确定度",
        "u(m/s)_标准不确定度",
        "u1(m/s)_瑞利探测通道光子噪声",
        "u2(m/s)_瑞利参考通道光子噪声",
        "u3(m/s)_瑞利探测通道背景噪声",
        "u4(m/s)_瑞利参考通道背景噪声",
        "u5(m/s)_瑞利探测通道渡越时间",
        "u6(m/s)_瑞利参考通道渡越时间",
        "u7(m/s)_瑞利测风锁频波长",
    ]
    result_df = pd.DataFrame(result, columns=header)
    # result_df.drop('U(%)_GUM扩展不确定度', axis=1, inplace=True)
    return result_df


def GUM_ResT(
    data_L2,
    data_L0,
    mode=0,
    c_T_RT=0.00349,
    c_T_RV=0.00349,
    table_T_resolution=0.1,
    tau=1e-9,
    u_tau=0,
    alt_angle=0,
    location_altitude=0,
    height_resolution=1,
    time_resolution=3600,
    pulse_repeat=30,
    k=1,
):
    # c_T_RT: 测温灵敏度，dT/dR_T
    # c_V_RV：测风灵敏度，dV/dR_V
    # data_L0：L0级数据，第1列为高度，第2列为RayVHS
    # data_L2：L2级数据，第1列为高度，第2列为密度TempRay
    # tau：探测器渡越时间
    # u_tau：探测器渡越时间波动度
    # pulse_repeat：激光器重频
    # time_resolution：L2级积分时间
    # location_altitude：系统所在高度
    # z0：参考密度所在高度
    # dens0：参考密度
    # uncertain_densityRef：参考密度波动度
    # tau = 1e-9
    # u_tau = 1e-9 * 0.01
    # time_resolution = 1200
    # pulse_repeat = 30
    # z0 = 30  #参考密度所在高度
    # dens0 = 3.64686762151e+017  #参考密度
    # u_dens0 = 0  #参考密度波动度
    # location_altitude = 1.14
    z = data_L2[:, 0]
    T = data_L2[:, 1]
    # 海拔高度转化和距离门内光子数累加
    data_L0_s = get_bins_data(
        data_L0, z, location_altitude, alt_angle, height_resolution
    )
    # 背景噪声获取
    B, B_mean_num = get_backgroud_noise(data_L0_s)
    B1 = B[0]
    B2 = B[1]
    B3 = B[2]
    # 探测范围截取
    Nr = get_target_data(data_L0_s, z)  # 目标高度光子数
    Nr1 = Nr[:, 0]
    Nr2 = Nr[:, 1]
    Nr3 = Nr[:, 2]
    # 光子数校正过程中的灵敏系数
    Nc1 = Nr1
    Nb1 = Nc1 - B1
    allSize = time_resolution * pulse_repeat  # L2数据累加的脉冲数
    c_Nb1_tau1 = Nc1**2 / allSize
    c_Nb1_Nr1 = Nc1**2 / Nr1**2
    c_Nb1_B = 1
    Nc2 = Nr2
    Nb2 = Nc2 - B2
    c_Nb2_tau2 = Nc2**2 / allSize
    c_Nb2_Nr2 = Nc2**2 / Nr2**2
    c_Nb2_B = 1
    Nc3 = Nr3
    Nb3 = Nc3 - B3
    c_Nb3_tau3 = Nc3**3 / allSize
    c_Nb3_Nr3 = Nc3**3 / Nr3**3
    c_Nb3_B = 1
    # 反演过程中灵敏系数
    if mode == 0:
        # 针对RT=(Nb1+Nb2)/Nb3/2, RV=(Nb1-Nb2)/Nb3
        c_RT_Nb1 = 1 / Nb3 / 2
        c_RV_Nb1 = 1 / Nb3
        c_RT_Nb2 = 1 / Nb3 / 2
        c_RV_Nb2 = -1 / Nb3
        c_RT_Nb3 = -(Nb1 + Nb2) / Nb3**2 / 2
        c_RV_Nb3 = -(Nb1 - Nb2) / Nb3**2
        # if c_T_RT == None:
        #     RT = (Nb1 + Nb2) / Nb3 / 2
        #     c_T_RT = np.polyfit(RT, T, 1)[0]
        # if c_T_RV == None:
        #     RV = (Nb1 - Nb2) / Nb3
        #     c_T_RV = np.polyfit(RV, T, 1)[0]
    c_T_RV = 0
    c_T_Nb1 = c_T_RT * c_RT_Nb1 + c_T_RV * c_RV_Nb1
    c_T_Nb2 = c_T_RT * c_RT_Nb2 + c_T_RV * c_RV_Nb2
    c_T_Nb3 = c_T_RT * c_RT_Nb3 + c_T_RV * c_RV_Nb3
    # 光子噪声不确定度分量
    u_Nr1 = np.sqrt(Nr1)
    u_Nr2 = np.sqrt(Nr2)
    u_Nr3 = np.sqrt(Nr3)
    c_Nr1 = c_T_Nb1 * c_Nb1_Nr1
    c_Nr2 = c_T_Nb2 * c_Nb2_Nr2
    c_Nr3 = c_T_Nb3 * c_Nb3_Nr3
    u1 = np.abs(u_Nr1 * c_Nr1)
    u2 = np.abs(u_Nr2 * c_Nr2)
    u3 = np.abs(u_Nr3 * c_Nr3)
    # 背景噪声不确定度分量
    u_B1 = np.sqrt(B1 / B_mean_num)
    u_B2 = np.sqrt(B2 / B_mean_num)
    u_B3 = np.sqrt(B3 / B_mean_num)
    c_B1 = c_T_Nb1 * c_Nb1_B
    c_B2 = c_T_Nb2 * c_Nb2_B
    c_B3 = c_T_Nb3 * c_Nb3_B
    u4 = np.abs(u_B1 * c_B1)
    u5 = np.abs(u_B2 * c_B2)
    u6 = np.abs(u_B3 * c_B3)
    # 渡越时间不确定度分量
    c_tau1 = c_T_Nb1 * c_Nb1_tau1
    c_tau2 = c_T_Nb2 * c_Nb2_tau2
    c_tau3 = c_T_Nb2 * c_Nb3_tau3
    u7 = np.abs(u_tau * c_tau1)
    u8 = np.abs(u_tau * c_tau2)
    u9 = np.abs(u_tau * c_tau3)
    # 温风查算表温度分辨率
    u10 = table_T_resolution / 2 / np.sqrt(3) * np.ones(z.shape)
    # 合成标准不确定度
    u = np.sqrt(
        u1**2
        + u2**2
        + u3**2
        + u4**2
        + u5**2
        + u6**2
        + u7**2
        + u8**2
        + u9**2
        + u10**2
    )
    # 扩展不确定度
    U = k * u
    # #不确定度
    U, u, u1, u2, u3, u4, u5, u6, u7, u8, u9, u10 = np.round(
        (U, u, u1, u2, u3, u4, u5, u6, u7, u8, u9, u10), 4
    )
    result = np.vstack((z, T, U, u, u1, u2, u3, u4, u5, u6, u7, u8, u9, u10)).T
    header = [
        "高度(km)",
        "温度(K)",
        "U(K)_GUM扩展不确定度",
        "u(K)_标准不确定度",
        "u7(K)_f0通道光子噪声",
        "u8(K)_f+通道光子噪声",
        "u9(K)_f-通道光子噪声",
        "u10(K)_f0通道背景噪声",
        "u11(K)_f+通道背景噪声",
        "u12(K)_f-通道背景噪声",
        "u13(K)_f0通道渡越时间",
        "u14(K)_f+通道渡越时间",
        "u15(K)_f-通道渡越时间",
        "u16(K)_温风查算表温度分辨率",
    ]
    result_df = pd.DataFrame(result, columns=header)
    return result_df


def GUM_ResW(
    data_L2,
    data_L0,
    mode=0,
    cRes_V_RT=0.00349,
    cRes_V_RV=0.00349,
    table_V_resolution=0.1,
    tau=1e-9,
    u_tau=0,
    alt_angle=0,
    location_altitude=0,
    height_resolution=1,
    time_resolution=3600,
    pulse_repeat=30,
    k=1,
):
    # c_T_RT: 测温灵敏度，dT/dR_T
    # c_V_RV：测风灵敏度，dV/dR_V
    # data_L0：L0级数据，第1列为高度，第2列为RayVHS
    # data_L2：L2级数据，第1列为高度，第2列为密度TempRay
    # tau：探测器渡越时间
    # u_tau：探测器渡越时间波动度
    # pulse_repeat：激光器重频
    # time_resolution：L2级积分时间
    # location_altitude：系统所在高度
    # z0：参考密度所在高度
    # dens0：参考密度
    # uncertain_densityRef：参考密度波动度
    # tau = 1e-9
    # u_tau = 1e-9 * 0.01
    # time_resolution = 1200
    # pulse_repeat = 30
    # z0 = 30  #参考密度所在高度
    # dens0 = 3.64686762151e+017  #参考密度
    # u_dens0 = 0  #参考密度波动度
    # location_altitude = 1.14
    z = data_L2[:, 0]
    V = data_L2[:, 1]
    # 海拔高度转化和距离门内光子数累加
    data_L0_s = get_bins_data(
        data_L0, z, location_altitude, alt_angle, height_resolution
    )
    # 背景噪声获取
    B, B_mean_num = get_backgroud_noise(data_L0_s)
    B1 = B[0]
    B2 = B[1]
    B3 = B[2]
    # 探测范围截取
    Nr = get_target_data(data_L0_s, z)  # 目标高度光子数
    Nr1 = Nr[:, 0]
    Nr2 = Nr[:, 1]
    Nr3 = Nr[:, 2]
    # 光子数校正过程中的灵敏系数
    Nc1 = Nr1
    Nb1 = Nc1 - B1
    allSize = time_resolution * pulse_repeat  # L2数据累加的脉冲数
    c_Nb1_tau1 = Nc1**2 / allSize
    c_Nb1_Nr1 = Nc1**2 / Nr1**2
    c_Nb1_B = 1
    Nc2 = Nr2
    Nb2 = Nc2 - B2
    c_Nb2_tau2 = Nc2**2 / allSize
    c_Nb2_Nr2 = Nc2**2 / Nr2**2
    c_Nb2_B = 1
    Nc3 = Nr3
    Nb3 = Nc3 - B3
    c_Nb3_tau3 = Nc3**3 / allSize
    c_Nb3_Nr3 = Nc3**3 / Nr3**3
    c_Nb3_B = 1
    # 反演过程中灵敏系数
    if mode == 0:
        # 针对RT=(Nb1+Nb2)/Nb3/2, RV=(Nb1-Nb2)/Nb3
        c_RT_Nb1 = 1 / Nb3 / 2
        c_RV_Nb1 = 1 / Nb3
        c_RT_Nb2 = 1 / Nb3 / 2
        c_RV_Nb2 = -1 / Nb3
        c_RT_Nb3 = -(Nb1 + Nb2) / Nb3**2 / 2
        c_RV_Nb3 = -(Nb1 - Nb2) / Nb3**2
        # if c_T_RT == None:
        #     RT = (Nb1 + Nb2) / Nb3 / 2
        #     c_T_RT = np.polyfit(RT, T, 1)[0]
        # if c_T_RV == None:
        #     RV = (Nb1 - Nb2) / Nb3
        #     c_T_RV = np.polyfit(RV, T, 1)[0]
    c_T_Nb1 = cRes_V_RT * c_RT_Nb1 + cRes_V_RV * c_RV_Nb1
    c_T_Nb2 = cRes_V_RT * c_RT_Nb2 + cRes_V_RV * c_RV_Nb2
    c_T_Nb3 = cRes_V_RT * c_RT_Nb3 + cRes_V_RV * c_RV_Nb3
    # 光子噪声不确定度分量
    u_Nr1 = np.sqrt(Nr1)
    u_Nr2 = np.sqrt(Nr2)
    u_Nr3 = np.sqrt(Nr3)
    c_Nr1 = c_T_Nb1 * c_Nb1_Nr1
    c_Nr2 = c_T_Nb2 * c_Nb2_Nr2
    c_Nr3 = c_T_Nb3 * c_Nb3_Nr3
    u1 = np.abs(u_Nr1 * c_Nr1)
    u2 = np.abs(u_Nr2 * c_Nr2)
    u3 = np.abs(u_Nr3 * c_Nr3)
    # 背景噪声不确定度分量
    u_B1 = np.sqrt(B1 / B_mean_num)
    u_B2 = np.sqrt(B2 / B_mean_num)
    u_B3 = np.sqrt(B3 / B_mean_num)
    c_B1 = c_T_Nb1 * c_Nb1_B
    c_B2 = c_T_Nb2 * c_Nb2_B
    c_B3 = c_T_Nb3 * c_Nb3_B
    u4 = np.abs(u_B1 * c_B1)
    u5 = np.abs(u_B2 * c_B2)
    u6 = np.abs(u_B3 * c_B3)
    # 渡越时间不确定度分量
    c_tau1 = c_T_Nb1 * c_Nb1_tau1
    c_tau2 = c_T_Nb2 * c_Nb2_tau2
    c_tau3 = c_T_Nb2 * c_Nb3_tau3
    u7 = np.abs(u_tau * c_tau1)
    u8 = np.abs(u_tau * c_tau2)
    u9 = np.abs(u_tau * c_tau3)
    # 温风查算表温度分辨率
    u10 = table_V_resolution / 2 / np.sqrt(3) * np.ones(z.shape)
    # 合成标准不确定度
    u = np.sqrt(
        u1**2
        + u2**2
        + u3**2
        + u4**2
        + u5**2
        + u6**2
        + u7**2
        + u8**2
        + u9**2
        + u10**2
    )
    # 扩展不确定度
    U = k * u
    # #不确定度
    U, u, u1, u2, u3, u4, u5, u6, u7, u8, u9, u10 = np.round(
        (U, u, u1, u2, u3, u4, u5, u6, u7, u8, u9, u10), 4
    )
    result = np.vstack((z, V, U, u, u1, u2, u3, u4, u5, u6, u7, u8, u9, u10)).T
    header = [
        "高度(km)",
        "风速(m/s)",
        "U(m/s)_GUM扩展不确定度",
        "u(m/s)_标准不确定度",
        "u8(m/s)_f0通道光子噪声",
        "u9(m/s)_f+通道光子噪声",
        "u10(m/s)_f-通道光子噪声",
        "u11(m/s)_f0通道背景噪声",
        "u12(m/s)_f+通道背景噪声",
        "u13(m/s)_f-通道背景噪声",
        "u14(m/s)_f0通道渡越时间",
        "u15(m/s)_f+通道渡越时间",
        "u16(m/s)_f-通道渡越时间",
        "u17(m/s)_温风查算表风速分辨率",
    ]
    result_df = pd.DataFrame(result, columns=header)
    return result_df


def GUM_RamanT(
    data_L2,
    data_L0,
    poly_coef,
    mode=0,
    tau=1e-9,
    u_tau=0,
    alt_angle=0,
    location_altitude=0,
    height_resolution=1,
    time_resolution=3600,
    pulse_repeat=30,
    k=1,
):
    # 输入量
    poly_coef = [0.003945, -0.00171, -0.00101]
    a = poly_coef[0]
    b = poly_coef[1]
    c = poly_coef[2]
    z = data_L2[:, 0]
    T = data_L2[:, 1]
    # 海拔高度转化和距离门内光子数累加
    data_L0_s = np.copy(data_L0)
    data_L0_s[:, 0] = (
        data_L0_s[:, 0] * np.cos(alt_angle / 180 * np.pi) + location_altitude
    )
    dz = np.abs(np.diff(z).mean())
    data_L0_s[:, 0] = np.round((data_L0_s[:, 0] - z[0]) / dz) * dz + z[0]
    data_L0_s = pd.DataFrame(data_L0_s).groupby(0, as_index=False).sum().values
    # 背景噪声获取
    B1 = data_L0_s[-1, 1]
    B2 = data_L0_s[-1, 2]
    # 探测范围截取
    ind = np.nanargmin(np.abs(z - z0))
    Nr2 = data_L0_s[ind, 1]  # 参考高度处光子数
    mask = (z[0] <= data_L0_s[:, 0]) & (data_L0_s[:, 0] <= z[-1])
    data_L0_s = data_L0_s[mask, :]
    Nr1 = data_L0_s[:, 1]
    # 光子数校正过程中的灵敏系数
    Nc1 = Nr1
    Nb1 = Nc1 - B1
    allSize = time_resolution * pulse_repeat  # L2数据累加的脉冲数
    c_Nb1_tau1 = Nc1**2 / allSize
    c_Nb1_Nr1 = Nc1**2 / Nr1**2
    c_Nb1_B1 = 1
    Nc2 = Nr2
    Nb2 = Nc2 - B2
    c_Nb2_tau2 = Nc2**2 / allSize
    c_Nb2_Nr2 = Nc2**2 / Nr2**2
    c_Nb2_B2 = 1
    # 反演过程中灵敏系数
    Q = Nb1 / Nb2
    temp = 1 / (c * np.log(Q) ** 2 + b * np.log(Q) + a)
    c_T_Nb1 = 1 / Nc2 * temp**2 * (2 * c * np.log(Q) + b) / Q
    c_T_Nb2 = Nc1 / Nc2**2 * temp**2 * (2 * c * np.log(Q) + b) / Q
    # 光子噪声不确定度分量
    u_Nr1 = np.sqrt(Nr1)
    u_Nr2 = np.sqrt(Nr2)
    c_Nr1 = c_T_Nb1 * c_Nb1_Nr1
    c_Nr2 = c_T_Nb2 * c_Nb2_Nr2
    u1 = np.abs(u_Nr1 * c_Nr1)
    u2 = np.abs(u_Nr2 * c_Nr2)
    # 背景噪声不确定度分量
    u_B1 = np.sqrt(B1)
    u_B2 = np.sqrt(B2)
    c_B1 = c_T_Nb1 * c_Nb1_B1
    c_B2 = c_T_Nb2 * c_Nb2_B2
    u3 = np.abs(u_B1 * c_B1)
    u4 = np.abs(u_B2 * c_B2)
    # 渡越时间不确定度分量
    c_tau1 = c_T_Nb1 * c_Nb1_tau1
    c_tau2 = c_T_Nb2 * c_Nb2_tau2
    u5 = np.abs(u_tau * c_tau1)
    u6 = np.abs(u_tau * c_tau2)
    # 合成标准不确定度
    u = np.sqrt(u1**2 + u2**2 + u3**2 + u4**2 + u5**2 + u6**2)


def get_bins_data(data_L0, z, location_altitude=0, alt_angle=0, height_resolution=1):
    # z = np.round(z, 3)
    data_L0_s = np.copy(data_L0)
    # plt.close()
    # plt.figure()
    # plt.semilogx(data_L0[:, 1:], data_L0[:, 0])
    # plt.show()
    data_L0_s[:, 0] = (
        data_L0_s[:, 0] * np.sin(alt_angle / 180 * np.pi) + location_altitude
    )  # 探测距离传化为海拔高度
    # dz = np.abs(np.diff(z).mean())
    # dz_L0 = np.abs(np.diff(data_L0_s[:, 0]).mean())
    dz_L0 = np.abs(data_L0_s[1, 0] - data_L0_s[0, 0])
    points = np.round(height_resolution / dz_L0)
    df_L0 = pd.DataFrame(data_L0_s)
    end_point = int(len(df_L0.index) // points * points)
    data_L0_s = (
        df_L0.iloc[:end_point, :]
        .groupby(np.arange(end_point) // points, as_index=False)
        .sum()
        .values
    )
    data_L0_s[:, 0] = data_L0_s[:, 0] / points

    # data_L0_s[:, 0] = np.round((data_L0_s[:, 0] - z[0]) / dz) * dz + z[0]
    # data_L0_s[:, 0] = np.round(data_L0_s[:, 0], 3)
    # data_L0_s = pd.DataFrame(data_L0_s).groupby(0, as_index=False).sum().values
    # data_L0_s[:, 0] = np.round(data_L0_s[:, 0], 3)
    # if height_resolution is not None:
    #     data_L0_s[:, 1:] = data_L0_s[:, 1:] * height_resolution / dz
    # plt.close()
    # plt.figure()
    # plt.semilogx(data_L0_s[:, 1:], data_L0_s[:, 0])
    # plt.show()
    return data_L0_s


def get_backgroud_noise(data_L0_s):
    mask = (130 <= data_L0_s[:, 0]) & (data_L0_s[:, 0] <= 160)
    B_mean_num = np.size(mask)
    B = np.nanmean(data_L0_s[mask, 1:], axis=0)  # 背景噪声取法1,以130km到160km均值为背景噪声
    # B = data_L0_s[-1, 1]  # 背景噪声取法2,以最后一个距离门为背景噪声
    return B, B_mean_num


def get_target_data(data_L0_s, z):
    # f = interpolate.interp1d(data_L0_s[:, 0], data_L0_s[:, 1:].T)
    # Nr = f(z).T
    f = interpolate.interp1d(data_L0_s[:, 0], np.log(data_L0_s[:, 1:].T))
    Nr = np.exp(f(z).T)
    return Nr
